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Abstract 

We present a new polarizable force field for aqueous ions (Li + , Na + , K + , Rb + , Cs + , Mg 2+ , 
Ca 2+ , Sr 2+ and CI - ) derived from condensed phase ab-initio calculations. We use Maximally 
Localized Wannier Functions together with a generalized force and dipole-matching procedure to 
determine the whole set of parameters. Experimental data is then used only for validation purposes 
and a good agreement is obtained for structural, dynamic and thermodynamic properties. The 
same procedure applied to crystalline phases allows to parametrize the interaction between cations 
and the chloride anion. Finally, we illustrate the good transferability of the force field to other 
thermodynamic conditions by investigating concentrated solutions. 



I. INTRODUCTION 



The development of classical force fields for ions in aqueous solution is essential to the 
description of specific effects, which are legion in biochemistry^-, atmospheric chemistry^ 
or environmental science^. The reliability of molecular simulations strongly depends on the 
quality of the force field used to represent the interactions, which must capture not only 
the effects of ionic size, but also the polarization of water by the ionic charge. The latter 
multi-body effect becomes essential when dealing with multivalent ions^, in concentrated 
solutions^ and in interfacial environments^^. 

In recent years, a successful strategy to derive polarizable force fields for solid and molten 
oxides from ab-initio simulations has been developed by Madden and co-workers^— . A full 
set of parameters was obtained for the Ca-Mg-Al-Si-0 (CM AS) system^, which is the main 
component of the Earth's crust and mantle. Cation-rich aluminosilicates, including clays 
and zeolites, are also the principal minerals on the Earth's surface, where they are in contact 
with ionic solutions. Examples of situations where the interface between such minerals and 
solution play an important role include the crystallization and dissolution of ionic crystals, 
such as calcium carbonate (in the context of carbon dioxide sequestration^) or sodium 
sulfate (deterioration of monuments^), or the sorption of radioactive contaminants (e.g. 
cesium or strontium) onto clays^ 1 ^. It is therefore of primary importance to extend the 
CMAS force field in order to describe these minerals and their interaction with water and 
ions. As a first step in this direction, we develop here a polarizable force field for ions in 
water which is compatible with that developed for the CMAS system. 

As mentioned previously, and despite the success of non-polarizable water force fields 
in the bulk— ~— , transferability to interfaces, especially charged ones, requires resorting to 
a polarizable model. Many such models exist, which differ mainly in their treatment of 
the polarizability. Drude or shell models assign a charge on a spring to each polarizable 
atom22~— . Other approaches allow for charge fluctuations^ 1 ^ or assign point dipoles to 
each polarizable species^— . Only the latter model is compatible with the above-mentioned 
one for oxides. Among the point polarizability models, we chose the one of Dang and Chang 
which was specifically developed to describe the gas-liquid interface^. In addition, Masia 
et al. have shown that it accurately reproduces the strong water polarization by divalent 
cations^^. 
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Following the strategy of Madden and co-workers, which has proven able to simultane- 
ously reproduce structural, dynamic and thermodynamic properties not only for the CMAS 
system, but also for many other ionic materials^^^i'^, we derive here the parameters of a 
force field for the aqueous ions: (Li+ Na+, K+ Rb+, Cs+, Mg 2+ , Ca 2+ , Sr 2+ and CI". The 
route from condensed phase Density Functional Theory (DFT) calculations, using Max- 
imally Localized Wannier Functions (MLWFs)^^ together with a generalized force and 
dipole-matching procedure^^, renders experimental input unnecessary, contrary to many 
force field parametrizations. 

The paper is organized as follows: We first give a detailed description of the force field and 
its parametrization which involves DFT calculations on single ions in bulk water and on ionic 
crystals. The second part is then devoted to the validation of the model, against structural, 
dynamic and thermodynamic properties of these systems. Finally, the transferability of the 
model is illustrated by the study of concentrated salt solutions. 

II. THE FORCE FIELD AND ITS PARAMETRIZATION 
A. Model 

The total energy of the system is decomposed into four terms: 

Vtot = Charge + Klisp + V rcp + V po \ (1) 

For the calculation of the direct Coulomb interaction between two atoms I and J, 

^charge = ^ 1 (2) 
I,J>I ' 

formal charges (here —1, +1 or +2) are used. The dispersion potential includes the dipole- 
dipole and dipole-quadrupole terms 



i,j>i 



/ 6 /J (^)T + /s J (^) 



r IJ r u 



(3) 



and the short-range corrections are described using the Tang-Toennies functions f™, which 
are of the form^: 

fl J = l-e^'»±®fe£ (4) 



k\ 

k=0 
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While the repulsion potential is modelled using a decaying exponential: 



[IJ e -B IJ ru 



(5) 



i,j>i 



Finally, many-body electrostatic effects are described by the induced dipoles fi 1 , which are 
treated as additional degrees of freedom and obtained at each MD step by minimizing the 
polarization energy: 

^poi = J2 ^7 W\ + E [W*W J {ru) - q'^iru)) Tfj ~ M T u] (6) 

with a 1 the ion polarizability and where the Einstein summation convention is assumed. A 
short-range correction to the multipolar expansion of the Tang-Toennies type is used: 



1 _ (JJe-^n 



k=0 



(b IJ r u y 
k\ 



(7) 



This so-called Polarizable Ion Model (PIM) has proven extremely successful for the 
description of oxides, chloride and fluoride-based materials, both in the solid and liquid 
statesii^ 1 ^. Water is described by a model compatible with this form, developed by Dang 
and Chang 2 ^. The only differences with the PIM are the description of the repulsive and dis- 
persion terms V rcp + Vdi sp for the water-water interactions, represented by a Lennard- Jones 
potential, and the absence of short-range damping of the charge-dipole interaction. The 
Dang-Chang (DC) water is a rigid 4-site model, with an additional virtual site M along the 
symmetry axis of the molecule, which bears a negative partial charge, as well as the induced 
dipole, while the Lennard- Jones interaction acts on the oxygen atom only. The parameters 
of the DC model are summarized in Table [H 



d H A 


doM A 


angle (°) 


eo (kcal/mol) 


a Q A 


QH 


«M A 


0.9752 


0.215 


104.52 


0.1825 


3.2340 


0.5190 


1.444 



TABLE I: Parameters of the Dang-Chang water model. 



The purpose of the present work is to derive all the parameters of the PIM for water-ion 
and ion-ion interactions, thereby providing a force field for the simulation of ions which is 
transferable from infinite dilution to concentrated solutions, up to the ionic crystals, for 
alkaline (Li+, Na + , K+, Rb + , Cs + ) and alkaline earth (Mg 2+ , Ca 2+ , Sr 2+ ) cations and the 

4 



chloride (CI - ) anion. Overall, this requires specifying 241 parameters. The procedure to 
determine all of them from ab-initio calculations aims at minimizing the risk of compensation 
of errors among the different terms by 1) directly computing as many parameters as possible, 
2) adjusting the remaining ones on different quantities (dipoles and forces) and 3) resorting 
to simplifying assumptions when necessary. We now describe these three aspects. 

B. Calculating parameters 

First-principle calculations based on Density Functional Theory (DFT) describe the elec- 
tronic density using the Kohn-Sham orbitals, whose delocalized nature renders the assign- 
ment of atomic or molecular properties difficult. The concept of the maximally localized 
Wannier function (MLWF) provides a convenient framework to analyze atomic and molec- 
ular properties in the condensed phased. The Wannier functions are defined through a 
unitary transformation of the Kohn-Sham eigenvectors. MLWFs are contructed by choosing 
the phase so that it minimizes the spread of the Wannier function^. It was shown recently 
that MLWFs could be used to systematically derive both the polarizabilities a 1 and disper- 
sion parameters C\ 3 and C\ J of a PIM— >22. Figured] illustrates the electronic density around 
a Ca 2+ cation and two water molecules in bulk water, reconstructed from their respective 
Wannier orbitals. 

1. Polarizability 

In a closed shell system, each MLWF contributes two electrons, so that the atomic or 
molecular dipole can be computed (in atomic units) as 



where the sums run over atoms i belonging to fragment / and over MLWFs n G i whose 
center is localized in the vicinity of the nuclear position Ri, Zi is the charge of nucleus i and 
is the position of the center of the n-th MLWF. The polarizability may differ from that 
of the same ion in the gas phase because of environmental effects. It can be calculated by 
applying a small electric field £^ along each Cartesian direction a = x, y, z to the system, 
which induces dipole moments l^ 1 } Ie H N .. A convenient way to distinguish the effect of 




(8) 
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FIG. 1: Localized electronic density around a Ca 2+ cation and two water molecules in bulk water, 
reconstructed from their respective Wannier orbitals. The isodensity surfaces include 90 % and 
95 % of the corresponding densities, respectively. 

the applied field from that of the static fields caused by the permanent charge distributions 
of the molecules, is to think of the former as an optical field. Sfi 1 can then be seen as the 
net induced dipole oscillating at the optical frequency. The total field f 1 '^ on each atom is 

f,(a) = £ {a) + J2T IJ - 6fl J (R N ) (9) 

where T IJ is the dipole-dipole interaction tensor. The polarizability tensor of molecule / 
can then be obtained by inverting Eq. (Q : 

a 7 ^) = {F 1 )" 1 ■ U 1 (10) 

with the second-rank three dimensional tensors defined as : = fi'^ and 11^ = S^id 
More details about this approach can be found in Ref.— . 
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2. Dispersion: C% and Cs 



DFT calculations do not usually account for dispersion interactions, because the former 
describe the electronic ground state while the latter arise from correlated density fluctuations 
associated with excited states. The treatment of dispersion via non-local functionals has only 
recently been introduced, albeit at a high computational cost. Thus these interactions are 
generally added (if at all) as an a posteriori correction. Among the several methods that 
have been proposed for computing this correction, the method of Grimme^ and that of 
Silvestrelli^ seem to be the most popular. In this work we use the latter, which considers 
the dispersion interaction between all pairs of MLWF as follows. The long-range interaction 
between separated fragments of matter is calculated, following Andersson et al.— , as 



6e f j„ VPiO , i)p 2 (r 2 ) w 1 



E xc = a i a nq/o i/o / / dridr 2 ^= y — . x ^ (11) 



where piri) is the charge density of fragment i, m the electron mass and Vi the volume 
occupied by fragment i. For large separations R, this scales as E = —Cq/R 6 , where the 
Cq coefficient for the interaction between two MLFWs k and I can be computed as : 



Cf = ArxH dr,dr 2 (12) 

ri<r c 
r2<r' c 

The cut-off radius r c = (1.475 — 0.866 In 5)5 is chosen to correctly capture the limit of 
long-range perturbations in an electron gas^. The MLWFs, giving rise to densities p = w 2 , 
are assumed to be of the Slater form: 

o3/4 

wj\\r - rvll) = ^ e -^ s ^ r ~ r ^ (13) 



-(VS/Sn 



t3/2 



characterized solely by their spread S n = {w n \r 2 \w n ) — (w n \r\w n ) 2 and center r n . 

We have previously shown that the dispersion interaction between two ensembles of charge 
density fragments can be obtained from the averaged sum over pair interactions of MLWFs^. 
Assuming an isotropic distribution of MLWF centers around the nuclei / and J, at fixed 
distances, leads (to second leading-order) to Vdi sp = ~^2 n =6s^n J l r lji where the dispersion 
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coefficients are: 



6 



(14) 



fce/,ieJ 



(7 J 
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£ 5(4 + ^)^6 



(15) 



kai,ia,j 



where g^,; are the distances of the MLWF centers to their respective nuclei and Cg is 
computed for each pair of MLWFs according to Eq. [T2l The determination of the parameter 
b D in Eq. HJ for the short-range damping of the dispersion interaction, is detailed below. 

C. Dipole- and force- fitting 

Not all parameters of the force field can be derived systematically from the electronic den- 
sity. However, they can be determined numerically so as to best reproduce the atomic prop- 
erties calculted by DFT: the total dipoles (permanent plus induced) of ions and molecules 
and the forces acting on them. 

1. Damping of charge- dipole interaction 

The first step in our parametrization procedure is to determine the parameters involved 
in Eq. [7] for the short-range damping of the charge-dipole interaction. This is achieved by 
numerically adjusting these parameters so as to minimize the error on the dipoles calcu- 
lated using the classical force field, relative to the DFT ones on a number of representative 
configurations: 



Together with the polarizabilities, these parameters complete the description of the polar- 
ization potential V po \. 

2. Repulsion 

The parameters of the repulsive potential V iep in Eq. \5\ can then be obtained by a similar 
procedure as the one used for the dipoles, if the functional used for the DFT calculation 




(16) 
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does not include dispersion interactions (e.g. PBE or BLYP}^>^: 

-^classical jpDFT\\2 



xl 



1 ^classical jplJp 1 l\Z 

— E E — 1 1 rrDFT\ 12 

(17) 



iVconf iVatan, c 



By ajusting the parameters for the damping of the charge-dipole interaction and for the 
repulsion on different physical quantities (dipoles and forces, respectively), we limit the risk 
of having a compensation of errors between the different terms of the potential. 



D. Further considerations 



The water-ion interactions are parametrized by applying the procedure described above on 
configurations of a system containing a single ion in bulk water. For the ion-ion interactions, 
we use configurations of the experimentally stable crystal phase under normal conditions: 
NaCl structure for Li + , Na + , K + and Rb + , CsCl structure for Cs + , MgCl2 structure for 
Mg 2+ , and CaCl2 structure for Ca 2+ and Sr 2+ . The Cl-Cl interactions must be the same 
among the different crystals in order to ensure the consistency and transferability of our 
potentials. The parameters for the Cl-Cl repulsion are obtained for LiCl, in which they are 
the most prominent, and the corresponding values are then used for all crystals. The Cq 
and C$ parameters for the Cl-Cl dispersion interaction are obtained by averaging the values 
for the different crystals. 

For the cation-anion repulsion (see Eq. E}, the force-fitting procedure results in B param- 
eters that are very close to each other among the alkaline ions on the one hand, and among 
the alkaline earth ions on the other hand. For the sake of simplicity, we use only one value 
for this parameter for each ion series. The A parameters for the cation-anion repulsion are 
then readjusted to minimize Eq. [17] while keeping the B value fixed. The final values for 
A and the corresponding Xf were practically unchanged by this constraint, thus confirming 
the relevance of this choice. 

In order to further decrease the number of free parameter, the range of the short-range 
damping used for the cation-anion dispersion bp (see Eqs. [3] and [4]) is taken in most cases 
equal to that of the short-range repulsion B IJ . This assumption is not new^, and it is jus- 
tified by the notion that the long-range scaling of dispersion breaks down as the electronic 
fragments start overlapping, when the short-range repulsion comes into play. The damp- 
ing of the Cl-Cl dispersion is adjusted numerically so as to reproduce simultaneously the 
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experimental density of all crystals. For the largest cations, Cs + and Sr 2+ , a value slightly 
smaller than £? catlon - cl was needed to reproduce the experimental densities. Compared to 
the usual procedure of parametrizing a PIM from ab-initio simulations^ 1 ^, the systematic 
determination of the Cg and Cg coefficients and the assumption that bp = B dramatically 
reduce the number of parameters that need to be adjusted in order to reproduce the whole 
set of experimental densities. The damping parameter of the monovalent cation-water dis- 
persion interaction was chosen equal to that for the corresponding monovalent cation-Cl _ 
dispersion interaction 6^ n "° = b 1 ^ , since the water molecules and the Cl~ ions have ap- 
proximately the same size. As far as the divalent cations are concerned, the attractive force 
arising at short distances from dispersion is negligible compared to the charge-charge and 
charge-dipole interactions. We can thus omit damping this interaction without any loss of 
accuracy. Similarly, the dispersion interaction between Cl~ and water oxygen atom is not 
damped. 

Overall, these considerations reduce the number of parameters for the interaction of all 
ions with water and of cations with chloride from 241 to 187, after the neglect of some terms 
for the reasons explained above, and to 170 by further assuming that the ranges of some 
interactions are equal. Out these 170, only 82 are adjusted numerically from the dipole- 
and force-matching procedures of section III CI while the rest are computed as expained in 
section III Bl 

E. Simulation details 

The parametrization of the force field from ab-initio simulations is achieved using repre- 
sentative configurations of the aqueous ions and the ionic crystals. For each ionic species, 
~ 100 configurations of a system containing a single ion and 32 water molecules are gener- 
ated using the force-field of Dang et al.—~— for the monovalent ions, and that of Yu et alM 
for the divalent ions. DFT calculations were then performed on these configurations with 
the BLYP functional^ 1 ^ (exept for the Rb + , for which the PBE functional^ was used). The 
Troullier-Martins^ (Cl~,Cs + and K+) and Goedecker-Teter-Hutter-^^ (Na+, Rb+, Mg 2+ , 
Ca 2+ and Sr 2+ ) pseudopotentials were used, with a plane-wave basis set and an energy cutoff 
of at least 70 Ry. Similarly, configurations of crystals containing between 16 and 108 MCI 
or MCI2 units, are used to perform the DFT calculations, with the same functionals and 
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pseudopotentials as for the ions in water. In each case, after determining the electronic 
density, the forces acting on each atom are computed and the dipoles are calculated from 
the MLWFs as described above. The C§ and C§ dispersion parameters are computed from 
the spreads and distances to the center of the MLWFs, which result from the localization 
procedure. The polarizabilities are calculated as explained above, by applying an external 
field using the Berry phase representation^ 2 -. All ab-initio calculations were performed using 
the CPMD simulation package^ (exept for those involving Ca 2+ , performed with CP2K 
simulation package^), while classical forces and dipoles are computed on the same con- 
figurations with FIST, the classical MD module of the CP2K simulation package^. The 
numerical minimization of Eqs. (TTB]) and ( II 7p is performed using the Minuit library^ 2 -. 

F. Parametrization: Results 

The computed polarizabilities for all the ions are summarized in table [III As expected, the 
polarizability increases when going down along columns of the periodic table (alkaline and 
alkaline earths series) , while a decrease is observed when going from left to right along rows 
(Na + to Mg 2+ , K + to Ca 2+ and Rb + to Sr 2+ ). For cations, the condensed phase polarizability 
is comparable to that in the gas phase, except for Cs + . For the chloride anion, however, the 
confinement of electrons by the surrounding water molecules results in a significant decrease 
of the polarizability (approximately 35%). A more detailed discussion has been given in 
RefA 

Interestingly, as indicated in Table [Til the (induced) dipole moment of cations is always 
very small compared to that of the chloride anion and water. This can be explained by the 
combination of two factors. First, most cations have a small (K + , Sr 2+ ) or very small (Li + , 
Na + , Mg 2+ , Ca 2+ ) polarizability. Second, all cations have a highly symmetric hydration 
sphere, which results in very weak local electric fields to polarize them. Because the induced 
dipoles are very small, they are not easily reproduced by the classical force field (typical 
errors are of the order of 100%), but they do not contribute significantly to the polarization 
energy V po i, which is dominated by the interaction of the ionic charge with the dipole of 
water, and hence to the forces. For the sake of simplicity, we thus decided to neglect the 
polarizability of all cations and not include any additional degrees of freedom to describe 
their induced dipoles. 
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Ion 


Q (A 3 ) 


V(/^ 2 ) (Debye) 


Li+ 


0.03 


0.002 


Na+ 


0.18 


0.014 


K+ 


0.81 


0.062 


Rb+ 


1.32 


0.097 


Cs+ 


2.02 


0.153 


Mg 2+ 


0.08 


0.010 


Ca 2 + 


0.44 


0.026 


Sr 2 + 


0.81 


0.071 


cr 


3.50 


0.415 



TABLE II: Polarizability a and magnitude of the induced dipole of each ion w (fJ> 2 )- The latter is 
1.18 D for water. 



The parameters for the cation-water interaction are summarized in Table 1111} and those 
for the chloride-water interaction are given in Table IIVI Finally, all parameters for the ion- 
ion interactions are given in Table |V] The resulting repulsion potentials V rep between water 
and the various cations, plotted in Fig. [2j nicely reflect the expected increase in ionic size 
along the alkaline and alkaline earth series. Furthermore, a comparable repulsion is observed 
for isoelectronic species such as Na + and Mg 2+ , K + and Ca 2+ , and Rb + and Sr 2+ . In line 
with the polarizabilities, the Cq and Cg dispersion coefficients for the ion-water interaction 
increase along the alkaline and alkaline earths series, while a decrease is observed from left to 
right along rows of the periodic table. The same trends hold for the repulsion and dispersion 
interactions between the cations and the chloride anion. 

We now examine the performance of the force field in terms of reproducing the ab- 
initio dipoles and forces. Figure |3] illustrates the comparison between the forces on the 
ion calculated with the classical force field (without dispersion) and those obtained from 
the DFT calculations, for the Ca 2+ cation. From table IVIt the relative error of the force, 
yxf^, on the Ca 2+ ion, with respect to the DFT result, is approximately 23% . This can 
be considered as a good match, especially when comparing to the corresponding results 
obtained by using the Dang potential 3 - (with the same water model), which results in a 
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System 


A ion -° (Ha) 


gion-O (£-1) 


Cj? n -° (Ha.A 6 ) 


Cf n ~° (Ha.A 8 ) 


b# (A- 1 ) 


b ion-M (A -1 ) 


c ion-M 


Li + -water 


24.75 


4.094 


1.103xl0~ 2 


1.037xl0~ 2 


3.000 


4.011 


2.950 


Na + -water 


711.1 


5.061 


1.335X10" 1 


1.572X10" 1 


3.000 


1.562 


6.839X10- 1 


K + -water 


125.7 


3.735 


7.530X10" 1 


1.206 


3.000 


1.315 


4.623xl0 -1 


Rb + -water 


157.8 


3.656 


1.225 


2.267 


3.000 


1.248 


4.765xl0 -1 


Cs + -water 


269.4 


3.635 


2.040 


4.644 


1.800 


2.524 


2.948 


Mg 2+ -water 


65.67 


3.963 


6.408xl0~ 2 


7.23xl0~ 3 




3.963 


2.820 


Ca 2+ -water 


57.94 


3.327 


5.055X10" 1 


7.502xl0 -1 




3.327 


3.000 


Sr 2+ -water 


41.55 


2.991 


9.159X10" 1 


1.576 




2.991 


2.041 



TABLE III: Parameters for the cation-water interactions. As for water-water interactions, respul- 
sion and dispersion involve the oxygen atom, while electrostatic interactions involve the additional 
M site. The damping parameter bp for the dispersion interaction for the monovalent ions is cho- 
sen equal to that of the corresponding cation-chloride interaction (see text and table [V]). The 
electrostatic damping is between the water dipole and cation charge. 



System 


A ion -° (Ha) 


gion-O (A" 1 ) 


Cjf 11 - (Ha.A 6 ) 


Cjf n -° (Ha.A 8 ) 


b ion-H ( A -l) 


ion-H 


b ion-M (A -1) 


c ion-M 


Cl-water 


499.63 


3.560 


2.039 


4.296 


4.794 


1.093 


2.444 


-1.901 



TABLE IV: Parameters for the chloride-water interactions. The dipole damping is between the CI" 
and the water charges. For the reasons already explained, there is no damping of the dispersion. 



relative error of approximately 320%. Similarly, our results for Sr 2+ show a relative error 
of 36%, compared to 501% with the force field from the literature^, 49% vs. 66% for Na + , 
48% vs. 131% for Cs + and 53% vs. 104% for CI - . Overall, the forces on all ions are well 
reproduced by the present force field. The largest contributions to the relative error (see 
Eq. [T7|) correspond to the smaller forces. 

Table IVHI reports the x 2 values obtained on the crystals for the forces on both cations 
and anions, as well as for the dipole of the anions. Comparison with Table [VTl indicates that 
a similar accuracy is obtained for both the crystals and the ions in solution, suggesting that 
the force field should perform well under both conditions. This result is also encouraging 
from the point of view of the transferability and the possible prediction of the solubility of 
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System 


Ion pair IJ 


A IJ (Ha) 


P/ J (A" 1 ) 


Ci J (HaJ 


6 ) 


Ci J (Ha.A 8 ) 


b# (A" 1 ) 


b IJ (A- 1 ) 


c u 


LiCl 


Li+-Li+ 


481 Q 




6.958 


2.727x10" 


-4 


5.570xl0~ 10 


6.958 








Li+-cr 


±0.0\J 




3.000 


2.369x10" 


-2 


2.511xl0~ 2 


3.000 


3.128 


1.433 




cr-cr 


698.4 




3.777 


5.951 




12.85 


1.650 






NaCl 


Na + -Na+ 


1 701 x 1 (T 


-2 


4.965 


2.914x10" 


-2 


1.394xl0~ 2 


4.965 








Na+-Cr 


44 4^ 




3.000 


2.971x10" 


-1 


3. 785x10^ 


3.000 


2.775 


2.040 




cr-cr 


698.4 




3.777 


5.951 




12.85 


1.650 






KC1 


K+-K+ 


1 74 Q 




5.000 


7.172x10" 


-1 


9.260xl0 _1 


5.000 








K+-cr 






3.000 


1.973 




3.347 


3.000 


1.282 


9.059X10" 1 




cr-cr 


698.4 




3.777 


5.951 




12.85 


1.650 






RbCl 


Rb+-Rb+ 


1 235 x 1 (T 


-2 


3.485 


2.235 




3.908 


3.485 








Rb+-cr 


108.0 




3.000 


3.755 




7.223 


3.000 


1.460 


9.825X10" 1 




cr-cr 


698.4 




3.777 


5.951 




12.85 


1.650 






CsCl 


Cs+-Cs+ 


353.0 




3.782 


7.325 




18.64 


3.782 








Cs+-cr 


150.1 




3.000 


7.339 




16.96 


1.800 


1.541 


4.665X10" 1 




cr-cr 


698.4 




3.777 


5.951 




12.85 


1.650 






MgCl 2 


Mg+-Mg 2+ 


2 231 x 1 0" 


-1 


4.995 


1.095x10" 


-2 


4.066xl0~ 3 


4.995 








Mg 2+ -cr 






3.400 


1.471x10" 


-1 


2.102xl0 _1 


3.400 


2.886 


2.113 




cr-cr 


698.4 




3.777 


5.951 




12.85 


1.650 






CaCl 2 


Ca 2+ -Ca 2 + 


1 289x10" 


-1 


3.941 


3.274x10" 


-1 


3.456xl0 _1 


3.941 








Ca 2+ -Cr 


236 3 




3.400 


1.168 




1.883 


3.400 


2.052 


1.268 




cr-cr 


698.4 




3.777 


5.951 




12.85 


1.650 






SrCl 2 


Sr 2 +-Sr 2 + 


5.513 




4.735 


1.259 




1.867 


4.735 








Sr 2+ -cr 


269.4 




3.400 


2.697 




5.066 


2.400 


3.103 


2.939 




cr-cr 


698.4 




3.777 


5.951 




12.85 


1.650 







TABLE V: Parameters for the ion-ion interactions. 



these crystals. Comparison between tables IVHI and IVIIII illustrates the better performance 
of the present model compared to those of Dang and coworkers^""—. 

Neglecting the polarizability of cations does not prevent us from obtaining a good de- 



ll 




FIG. 2: Repulsion potential between water and alkaline cations (a) and alkaline earth cations (b), 
in units of the thermal energy /3" 1 = fcgT. 

scription of the forces acting on them, as can be seen in table I VII These forces are even 
better described than those on the chloride ion, whose polarizability is explicitly taken into 
account. Nevertheless, for the reasons mentioned in the introduction, it is essential to cor- 
rectly reproduce the polarization of water molecules around ions. Table I VI I also indicates 
the relative error on the dipole of water molecules in the first solvation shell of the ions. 
The combination of the Dang-Chang water model with the present model for the ion-water 
interactions provides a very good description of the polarization of water, with relative errors 
between 5 and 10% for all ions except Mg 2+ (13%). 

III. VALIDATION 

Having shown that our force field is able to correctly reproduce the ab-initio dipoles and 
forces, we now turn to its validation against experimental data pertaining to the structure, 
thermodynamics and dynamics of aqueous ions at infinite dilution, as well as to the density 
of ionic crystals. We finally investigate the transferability of the force field to concentrated 
solutions, which where not taken into account when "designing" the force field. It is worth 
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Configuration index 



FIG. 3: Force (in atomic units) acting on the Ca 2+ ion. The prediction of the classical force 
field (lines) for the 3 components F x , F y and F z are compared to the DFT result (x), for 100 
configurations composed of 32 water molecules and 1 Ca 2+ . 

pointing out here that we use experimental data only for validation purposes, in contrast 
with all other force fields for aqueous ions, which use some experimental data for calibration 
of the parameters. Out of the 241 parameters defining the force field for the present set 
of ions, only 3 (the dispersion damping parameters fo^ 1_cl , b^~ cl and 6§" C1 ) are determined 
with the use of experimental data, namely the densities of the 8 crystals. In particular, no 
experimental data on aqueous ions is used during the calibration process. 

A. Simulation details 

For ions at infinite dilution, the system contains a single ion and 215 water molecules in 
a cubic box of size L = 18.65 A. For the crystals, the systems consist of 256 LiCl, NaCl, KC1 
or RbCl, 342 CsCl, 192 MgCl 2 or CaCl 2 , or 256 SrCl 2 . Systems for concentrated solutions 
are composed of 27 NaCl, KC1 and 458 water molecules in cubic box of sizes 24.4167 A 
and 24.638 A, respectively. Electrostatic interactions are computed using a dipolar Ewald 
sum£L£&, with a tolerance of 1.10~ 7 to obtain the self-consistent dipole moments. Molecular 
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Ion 


^F-ion 


V 2 

A-lx - H 2 


X 2 ■ 

A,^£-ion 


Li+ 


1.56X10" 1 


1.76xl0" 3 




Na+ 


2.36X10" 1 


7.81 xlO" 3 




K+ 


l.llxlO" 1 


2.79x10-3 


_ 


Rb+ 


1.02xl0 _1 


2.35xl0" 3 


_ 


Cs+ 


2.28XKT 1 


3.08xl0" 3 


- 


Mg 2 + 


9.97xl0- 2 


1.78x10^2 


_ 


Ca 2 + 


5.35xl0- 2 


9.05X10" 3 




Sr 2 + 


1.27xl0 _1 


2.73X10" 3 




cr 


2.84X10" 1 


3.40 xlO" 3 


2.05xl0 _1 



TABLE VI: \ 2 f° r the forces on the ions and the dipoles of water and the ions. 



Crystal 


X F _ M x + 


X F-C1" 




LiCl 


1.13X10" 1 


1.74xl0- 2 


1.91X10- 1 


NaCl 


2.52x10^2 


1.12x10-2 


1.87X10- 1 


KC1 


7.88x10^2 


5.61xl0- 2 


7.28X10- 1 


RbCl 


4.77x10^2 


6.46x10-2 


6.66X10- 1 


CsCl 


2.23x10^2 


1.15X10- 1 


4.30X10- 1 


MgCl 2 


2.62xl0^ 1 


8.45x10^2 


1.90x10-2 


CaCl 2 


5.58x10^2 


8.46xl0- 3 


2.27X10- 1 


SrCl 2 


3.61x10^2 


3.95x10-2 


5.69x10-2 



TABLE VII: \ 2 hi crystals, for the forces on the cations and anions, and the dipoles of anions. 

dynamics in the canonical ensemble are performed using a Nose-Hoover thermostat with a 
time constant of 1 ps. The system is first equilibrated for 250 ps, and the properties are 
determined from subsequent 2.75 ns runs. The density of the crystals is determined from 
simulations in the NPT ensemble at P = 1 bar. The thermostat is the same than the one 
used for the NVT ensemble and the barostat is an extension of the one by Martyna et al.—. 
All simulations are performed using the CP2K simulation package^. 



17 



Crystal 


Y 2 + 


Y 2 

^F-Cl 


LiCl 


3.51 


28.0 


NaCl 


4.46XKT 1 


3.50 


KC1 


4.38 


4.10 


CsCl 


2.24 


2.56 


CaCl 2 


1.03 


9.85X10- 1 


SrCl 2 


5.47 


8.28 



TABLE VIII: y 2 m crystals, for the forces on the cations and anions, with the polarizable Dang- 
Chang models. 



B. Solvation of ions: structure 



We first inverstigate the structure of the solvation shells around ions by computing radial 
distribution functions, reported for the cations in Fig. HI As usually observed, the position of 
the first maximum gradually shifts towards larger distances when switching from Li + to Cs + 
and from Mg 2+ to Sr 2+ , while the value of the maximum decreases and the peak broadens. 
On the contrary, moving right along the rows of the periodic table results in a closer and 
sharper peak. This arises from the stronger electrostatic interaction, since the ion-water 
repulsion remains comparable, as discussed previously, and reflects a tighter first solvation 
shell. 




1 




i 1 i 

- M g 2+ ; 






- Ca 2+ 

sr 2+ ; 




v.. 











4 

r(A) 



FIG. 4: Ion-oxygen radial distribution functions for the aqueous cations. 



18 



The positions of the first maximum and the coordination numbers, defined as the integral 
of the ion-0 (water) radial distribution function from the origin out to the first minimum, 
are summarized in Table IIX[ together with the corresponding experimental values. The 
value and error estimates of the coordination numbers are determined from the plateau of 
the running values. Remarkably, all simulated data fall in the reported experimental ranges. 
Particularly encouraging is the agreement with experimental data for the three divalent ions. 
While several force fields are able to correctly predict the position and number of neighbours 
for the Mg 2+ ion, many of them fail to correctly reproduce that of Ca 2+ . As an example, a 
force field by Yu et al. based on a Drude model of polarizability, which accurately describes 
the hydration free energies, predicts a coordination number of 6 for this ion^ 2 -. Our result 
is very close to the value of 7.3 obtained with the popular AMOEBA force field used for 
biomolecular simulations^. The previously available model for Ca 2+ , with the present water 
model, predicts a distance of 2.45 A, within the experimental range, but it used the EXAFS 
data of 2.43 A in the parametrization process. In the case of Sr 2+ , we find a distance very 
close to the anomalous X-ray diffraction value of 2.67A— and a coordination number which 
is within the reported experimental range. 

When comparing simulation results for an ion at infinite dilution with experiments, one 
should pay attention to the experimental conditions, in particular the concentration and the 
nature of the counterion. For example, Smirnov and Trostin reported an increase of the 
Cs + coordination number with decreasing concentration^. It is thus not surprising to find 
our result on the larger side of the experimental range. Moreover, results using CIO4 as a 
counterion instead of Cl~ are less likely to be polluted by the formation of ion pairs. The 
distances of 2.12 and 2.65 A between the cation and the nearest water oxygen, reported with 
ClOj in Ref.— for Mg 2+ and Sr 2+ , respectively, are in very good agreement with ours (2.13 
and 2.68). 

Since most force fields include some experimental data on the structure during the cal- 
ibration process^ 2 - 1 ^ 1 ^, it is possible to obtain a good agreement. When such data is not 
included as a target property, the predicted structure may not be very accurate. As an 
example, Horinek et al. parametrized a simple non-polarizable force field optimized for the 
simulation of solvation thermodynamics^ 2 -. The structural properties, used only for vali- 
dation purposes, revealed a tendency to underestimate the distances to the nearest water 
molecules for cations. Our results for the chloride ion are very good, as they fall exactly on 
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the EXAFS value of 2.11 A determined by Dang et al, whereas many force fields tend to 
predict too large a distance for the first peak— i 7 - 2 -, even though they include such structural 
properties in the fitting procedure^ 3 -. 



Ion 


Position (A) 


Coordination Number 




Sim 


Exp 


Sim 


Exp 


Li+ 


1.96 


1.90-2.25 


4.0 


4 


Na+ 


2.41 


2.41-2.50 


5.7±0.1 


4-8 


K+ 


2.74 


2.60-2.92 


6.45±0.25 


4-8 


Rb+ 


2.88 


2.80-3.05 


7.05±0.25 


6-8 


Cs+ 


3.20 


2.95-3.21 


8.3±0.8 


6-8 


Mg 2+ 


2.13 


2.00-2.15 


6.0 


6 


Ca 2+ 


2.53 


2.40-2.58 


7.24±0.02 


7-9 


Sr 2 + 


2.68 


2.57-2.67 


7.81±0.05 


7.3-10.3 


cr 


3.11 


3.05-3.18 


6.12±0.12 


5.3-6.4 



TABLE IX: Structural properties: position of the first maximum in radial distribution function 
and coordination number. The experimental values are taken from Refs i 63 ' 69-71 ' 73 

Positions and coordination numbers cannot be measured directly, and the experimental 
values are the outcome of a complex numerical analysis of the raw data, which typically 
involves several Fourier transforms and filters which can influence the final result. A more 
stringent test of the force field thus consists in comparing the experimental signal to that ob- 
tained by computing the experimental observables on configurations generated by molecular 
simulation. An example of such a test is given in Fig. |5l which compares the experimental 
EXAFS signal for aqueous Ca 2+ , obtained from Ref.— , to that predicted from our configu- 
rations using the FEFF8 code, which uses an updated version of the Rehr et al. algorithm^ 
to evaluate multiple electron scattering series. The agreement is seen to be very good in the 
k > 3 A" 1 part of the spectrum, both in terms of the amplitude (which reflect the number 
of neighbours) and the frequency of oscillations (related to their position). 



20 



2 




5 



10 



k (A-i) 



FIG. 5: Comparison between simulated and experimental EXAFS^ signal for aqueous Ca 2+ . The 
agreement is very good both for the amplitude, which reflects the number of neighbours, and 
frequency of the oscillations, related to their position. 

C. Solvation of ions: hydration free energy 

Among all the ionic properties one aims to predict, the hydration free energy AGh y d 
is probably the most important, since it relates to the ability of ions to accomodate their 
solvation shell when approaching an interface or other ions. This quantity is almost always 
one of the target properties used to design force fields. Whereas absolute values are difficult 
to determine, differences in hydration free energies can be easily computed using a thermo- 
dynamic integration procedure without worrying about the numerous corrections^^ (for 
system size, boundary conditions, and the treatment of electrostatic interactions) needed for 
the former. The difference AAG hyd « AAF hyd = AF^ yd - AF^ A can be determined from 
a thermodynamic path (transmutation) connecting the systems, by introducing a mixed 
Hamiltonian H(X) = XH^ + (1 — X)H^ a for A G [0, 1], as 



For the monovalent ions, we use a 6-point Gaussian quadratures^ to compute the integral, 
except for the Li + -Cs + transmutation, for which we use an 8-point quadrature. Details 
on this standard quadrature procedure can be found in Ref.— . In the case of the divalent 
ions, where d\H(\) shows a linear variation in A, a simpler trapezoidal rule can be used to 




(18) 
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approximate the integral. In this case we used ten equally spaced points (0.1) for A, within 
the interval [0,1]. 



Transmutation 


AAGH (kcal/mol) 


AAGfP (kcal/mol) 


Li+^Na+ 

K+^Rb+ 
Rb+^Cs+ 


26.5 
13.7 
3.2 
7.6 


[23.8;26.2] 
riB 7-17 71 
14.9:5.41 
[5.5;7.7] 


Li+^Cs+ 


51.4 


[50.9;57.0] 


Mg 2 +^Ca 2 + 
Ca 2+ ^Sr 2 + 


82.2 
25.3 


[77.7;80.3] 
[29.8;32.9] 


Mg 2 +^Sr 2+ 


107.8 


[107.5;113.2] 



TABLE X: Differences in Gibbs free energy of hydration: Simulated and experimental values. The 
experimental values for the monovalent ions are taken from references' 8 -^ and those for divalent 
ions fro m 81 ' 83 ' . 



The hydration free energy differences, for all the transmutations considered, are summa- 
rized in Tab. |23> together with the corresponding experimental values. The overall agreement 
with experiment is very good, with deviations never exceeding a few kcal/mol, and the large 
variations of AAGh y d across the ion series being well reproduced. We note that some force 
fields are able to reproduce this quantity slightly more accurately, such as the non-polarizable 
one of Horinek et a/.— or the polarizable model (Drude oscillators) of Yu et alr^-. In these 
cases, however, experimental hydration free energies (or differences) were used as a target 
property to calibrate the force field, whereas we use it here as an independent validation of 
our ab-initio derived model. 



D. Diffusion coefficient 



The diffusion coefficient are computed using from the mean-squared displacement, as 
determined by the Einstein relation : 

l d(|r(t)-r(0)| 2 ) 

DpBC = ^6 at (19) 
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The "PBC" subscript emphasizes the fact that the use of periodic boundary conditions 
induces a box length dependence on the measured diffusion coefficient, which takes the 

f orm 86. 

DpB0 . B „_^™ (20 ) 



6nr]L 

where r] is the shear viscosity of the solvent. For the box length of L = 18.65 A used in our 
simulations, the correction to the Dang-Chang water model is approximately 0.43 10 -9 m 2 s _1 
and must not be neglected (D^ 2 ° = 2.72 ±0.09 m 2 s -1 )^. For a meaningful comparison with 
experiments, we thus extrapolate to the infinite box length limit, both for the ion and the 
water diffusion coefficients and compare the ratios D$ n /D$ 2 °. 



Ion 


(D$ n /D^ 2 °) sim 




Li+ 

Na+ 
K+ 
Rb+ 
Cs+ 


0.49 
0.54 
0.78 
0.88 
0.82 


0.44 
0.58 
0.88 
0.90 
0.89 


Mg 2 + 
Ca 2+ 
Sr 2 + 


0.31 
0.35 
0.35 


0.31 
0.34 
0.34 


cr 


0.71 


0.88 



TABLE XI: Ratio between the ion and water diffusion coefficients. The experimental values for 
the ions are taken from^, the one for the water from^ 9 -. 



The simulation results are compared to the experimental ratios in Tab. IXII The relative 
error is of only 2 to 11% for the monovalent cations. The agreement is particularly good for 
the divalent cations, for which the relative error does not exceed 3%. The largest relative 
error is for the Cl~ ion (19%). We performed similar simulations with the force field of Dang 
and co-workers^ 9 - - — , which uses the same water model. The errors in that case reach 31% 
for Cl~ and 9% for Ca 2+ and Sr 2+ . Moreover, while our results capture the equal diffusion 
coefficients of these two cations, the model of Dang and co-workers underestimates that of 
Ca 2+ and overestimates that of Sr 2+ . 
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Since force fields generally do not include experimental data on the dynamics for their 
calibration, their accuracy for dynamical properties is usually not as high as for structural 
and thermodynamic ones. The present strategy, which aims at reproducing the forces on the 
atoms and molecules, as best as possible, allows us to also predict the dynamic properties. 

E. Crystal density 

As a test of the interactions between ions, we now turn to the study of the ionic crystals. 
As explained above, out of the 241 parameters defining the force field for the entire family 
of ions we have studied, only 3 (the dispersion damping parameters fc^ 1 " 01 , 6§ S_C1 and fofj ) 
were determined using the experimental densities of the 8 crystals as target properties. All 
the systems we studied preserved their correct crystal structure during the entire length of 
the simulations, even the more complex ones corresponding to the divalent cations. Fig. O 
illustrates the deformed rutile structure of CaCl2 and the lamellar one of MgCl2. While a 
complete study of the relative stability of the different possible phases exceeds the scope of 
the present work, this suggests that these phases are at least metastable. The simulated 
densities are compared to the experimental ones in Tab. Kill The overall agreement is once 
again good, with relative errors below 10% except for NaCl (16%). 



Crystal 


p sim (g.cm -3 ) 


p exp (g.cm" 3 ) 


LiCl 


2.01 


2.07 


NaCl 


1.83 


2.17 


KC1 


1.93 


1.99 


RbCl 


2.98 


2.76 


CsCl 


4.42 


3.99 


MgCl 2 


2.21 


2.33 


CaCl 2 


2.04 


2.15 


SrCl 2 


3.25 


3.05 



TABLE XII: Density of the crystals at 1 bar and 300 K. The experimental values are taken 
£ rom 88^ Note that the correct crystal structures (separated in the table) are preserved during the 
simulations. 
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FIG. 6: A) Snapshot of CaCl2 cristal. B) Snapshot of MgCl2 cristal. Both structures are stable 
during the simulations. Cl~ are in cyan, Mg 2+ in red and Ca 2+ in green. 

Transferability to crystals is rarely tested, making comparisons with other potentials 
rather difficult. We have again used the force fields of Dang and coworkers^ - — to as- 
sess the reliability of our potentials. Although their potentials give good results for NaCl 
(d=2.1 g.cm -3 ), KC1 (1.9 g.cm -3 ) and CsCl (3.8 g.crn" 3 ) crystals, the structure proves to 
be unstable for LiCl, CaCl2 and SrCl2. This example shows the need for more complicated 
force fields (with more parameters), as they can provide better transferability. 

F. Concentrated solutions 

The previous sections demonstrate the accuracy of the present force field for both in- 
finitely dilute solutions and crystals. We now test its transferability to conditions which 
were not considered during the construction of the force field, by investigating concentrated 
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ionic solutions. We compute the neutron diffraction spectra for concentrated NaCl and KC1 
solutions with one ion pair for 17 water molecules (1:17) from the site-site partial structure 
factor between site a and (3: 

S a p(Q) = Anp j r\g aP {r) - V^^-dr (21) 

where p is the atomic number density of the solution and g a p{r) the corresponding site-site 
radial distribution function. Experimental neutron diffraction allows for the extraction of 
composite partial structure factors. We compare our simulations results to the traditionally 
used Fxx function, defined as: 

{Y, a c a b a f 

where the sums over a and (3 run over all atom types except hydrogen and c a and b a are the 
atomic fraction and neutron scattering length of atom a, respectively. The comparison with 
the experimental results taken from Ref.— in Figs. [7] and [8] indicates a very good agreement, 
which confirms the transferability to concentrated solutions. 



F XX (Q) = ^ V ^^jp^j (22) 



1 



§ 1 
-2 



7 
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o Experiment 



10 15 

Q (A 1 ) 



FIG. 7: Comparison between simulated and experimental Fxx(Q) from Ref.— for a concentrated 
NaCl solution (one NaCl pair for 17 water molecules). 



IV. CONCLUSION 

We have shown a successful parametrization of a polarizable force field for aqueous solu- 
tions of Li+, Na+, K+, Rb+, Cs+, Mg 2+ , Ca 2+ , Sr 2+ and CR ions. We used the polarizable 
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FIG. 8: Comparison between simulated and experimental Fxx(Q) from Ref.— for a concentrated 
KC1 solution (one NaCl pair for 17 water molecules). 

Dang-Chang model for water and derived all the parameters involving ions in the frame- 
work of the polarizable ion model of Madden and co-workers. The procedure relies only 
on ab-initio DFT calculations; part of the parameters (polarizabilities, dispersion coeffi- 
cients) are directly calculated while the others are extracted from a generalized force- and 
dipole-matching procedure. Experimental information is used for validation purposes only: 
The structural (first-neighbour distances, coordination numbers), thermodynamic (hydra- 
tion free energy differences) and dynamic (diffusion coefficients) are very well reproduced. 
The interactions between cations and the chloride anion are parametrized on calculations 
performed in the crystal phases, thus ensuring the accuracy of the force field across the 
whole concentration range. 

The account of multi-body effects via the polarizability should ensure a good transfer- 
ability to more complex conditions: mixtures of these salts, high temperature and pression^i 
and to liquid-vapor or liquid-solid interfaces. The next step will consist in extending the 
present approach to the interaction of water with the surface of oxide materials. It will 
then be possible to use this force field for the study of important problems of environmental 
science such as the retention of radionuclides onto clay minerals, or the water uptake by 
clays and zeolites. 
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